home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / procssng / ccs / ccs-11tl.lha / lbl / hips / sources / vfft / inv-vfft.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-03-12  |  8.0 KB  |  344 lines
  1. /*
  2. %    INV-VFFT . C
  3. %
  4. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  5.  
  6. This software is copyright (C) by the Lawrence Berkeley Laboratory.
  7. Permission is granted to reproduce this software for non-commercial
  8. purposes provided that this notice is left intact.
  9.  
  10. It is acknowledged that the U.S. Government has rights to this software
  11. under Contract DE-AC03-765F00098 between the U.S.  Department of Energy
  12. and the University of California.
  13.  
  14. This software is provided as a professional and academic contribution
  15. for joint exchange. Thus, it is experimental, and is provided ``as is'',
  16. with no warranties of any kind whatsoever, no support, no promise of
  17. updates, or printed documentation. By using this software, you
  18. acknowledge that the Lawrence Berkeley Laboratory and Regents of the
  19. University of California shall have no liability with respect to the
  20. infringement of other copyrights by any part of this software.
  21.  
  22. For further information about this notice, contact William Johnston,
  23. Bld. 50B, Rm. 2239, Lawrence Berkeley Laboratory, Berkeley, CA, 94720.
  24. (wejohnston@lbl.gov)
  25.  
  26. For further information about this software, contact:
  27.     Jin Guojun
  28.     Bld. 50B, Rm. 2275, Lawrence Berkeley Laboratory, Berkeley, CA, 94720.
  29.     g_jin@lbl.gov
  30.  
  31. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  32. %
  33. %    Inverse Virtual Fast Fourier Transform
  34. %
  35. %    Only one write_header() between header processing and inv_vfft.
  36. %
  37. % AUTHOR:    Jin Guojun - LBL    5/10/91
  38. */
  39. char    usage[]="options\n\
  40. -B    output BYTE format. Default is float.\n\
  41. -D    2D test    \n\
  42. -cut    cut both ends and output BYTE image\n\
  43. [<] VFFT [> Image]\n";
  44.  
  45. #include "complex.h"
  46. #include "header.def"
  47. #include "imagedef.h"
  48.  
  49. U_IMAGE    uimg;
  50.  
  51. bool    cut;
  52. MType    dimen1len, dimen2size, fsize, vsize;
  53.  
  54. #define    row    uimg.height
  55. #define    cln    uimg.width
  56. #define    frm    uimg.frames
  57.  
  58.  
  59. main(argc, argv)
  60. int    argc;
  61. char*    argv[];
  62. {
  63. bool    dimens=0, vflag;
  64. int    i, j, f, hrows, hcols, logfrms, logrows, logcols;
  65.  
  66. format_init(&uimg, IMAGE_INIT_TYPE, HIPS, -1, *argv, "S20-1");
  67. uimg.o_form = IFMT_FLOAT;
  68.  
  69. for (i=1;i<argc;i++) {
  70.     if (argv[i][0] == '-'){
  71.     switch(argv[i][1]) {
  72.     case 'c':    cut++;
  73.     case 'B':    uimg.o_form = IFMT_BYTE;
  74.         break;
  75.     case 'D':    dimens++;
  76.         break;
  77.     default:
  78. info:        usage_n_options(usage, i, argv[i]);
  79.     }
  80.     }
  81.     else if (freopen(argv[i], "r", stdin) == NULL)
  82.     syserr("can't open %s as input", argv[i]);
  83. }
  84.  
  85. io_test(stdin_fd, goto    info);
  86.  
  87. (*uimg.header_handle)(HEADER_READ, &uimg, 0, 0);
  88.  
  89. vflag = uimg.in_form==IFMT_DVVFFT3D || uimg.in_form==IFMT_VVFFT3D;
  90. if (uimg.o_form==IFMT_BYTE)
  91.     uimg.pxl_out = 1;
  92. else    uimg.pxl_out = 4;
  93.  
  94. hrows = row >> 1;
  95. hcols = cln >> 1;
  96. dimen1len = hcols+1;
  97. vsize = row * dimen1len;
  98. fsize = row * cln;
  99. if (!vflag)    /* vvfft must be in 3D. */
  100.     dimens |= !(uimg.in_form & 1);
  101.  
  102. logfrms=logrows=logcols = -1;
  103. for (i=0,j=1;i<12;i++,j+=j) {
  104.     if (j==row)    logrows=i;
  105.     if (j==cln)    logcols=i;
  106.     if (j==frm)    logfrms=i;
  107. }
  108.  
  109. (*uimg.header_handle)(HEADER_WRITE, &uimg, argc, argv, True);
  110.  
  111. if (logfrms<0 || logrows<0 || logcols<0){
  112.     mesg("not power of 2, be patient for slow processing\n");
  113.  
  114. #    include    "inrlvfft.cxx"
  115.     exit(0);
  116. }
  117.  
  118. w_init(MAX(MAX(hrows, hcols), 1<<logfrms-1));
  119.  
  120. if (uimg.in_form==IFMT_DVFFT3D || uimg.in_form==IFMT_DVFFT2D){
  121. double    *temp, *re_i, *im_i, *cvt,
  122.     *re_o = (double*)nzalloc(fsize, sizeof(double)<<1, "ibuf"),
  123.     *im_o = re_o + fsize;
  124.  
  125.     if (!dimens){
  126.         re_i = (double*)nzalloc(frm*vsize, sizeof(double)<<1, "obuf");
  127.         im_i = re_i + frm*vsize;
  128.         cvt = re_o;
  129.     }
  130.     else    cvt = (double*)nzalloc(fsize, sizeof(*cvt)<<1, "cvt");
  131.  
  132.     if (!dimens){
  133.         w_load(1<<logfrms);
  134.         if (vflag){
  135.             for (f=0; f<frm; f++){
  136.             upread(re_i+f*vsize, sizeof(*re_i), vsize, stdin);
  137.             upread(im_i+f*vsize, sizeof(*im_i), vsize, stdin);
  138.             }
  139.             for (i=frm*vsize; i--;)
  140.             im_i[i] = -im_i[i];
  141.         }
  142.         else    {
  143.         register double    *fop,
  144.             *re_c = re_i,
  145.             *im_c = im_i;
  146.             for (f=0; f<frm; f++){
  147.                 fop = cvt;
  148.                 i = upread(fop, sizeof(*fop)<<1, vsize, stdin);
  149.                 if (i != vsize)
  150.                 syserr("f%d [%d] read %d", f, vsize, i);
  151.                 for (i=vsize; i--;){
  152.                 *re_c++ = *fop++;
  153.                 *im_c++ = -*fop++;
  154.                 }
  155.             }
  156.         }
  157.         for (i=0; i<vsize; i++)
  158.             dvfftn(re_i+i, im_i+i, logfrms, vsize);
  159.     }
  160.  
  161.     for (f=0; f<frm; f++){
  162.         if (dimens){
  163.         register double    *re_out = re_o, *im_out = im_o;
  164.             temp = cvt;
  165.             upread(temp, sizeof(*temp)<<1, vsize, stdin);
  166.             for (i=0; i<row; i++){
  167.                 for (j=0; j<dimen1len; j++){
  168.                 *re_out++ = *temp++;
  169.                 *im_out++ = -*temp++;
  170.                 }
  171.                 re_out += cln - dimen1len;
  172.                 im_out += cln - dimen1len;
  173.             }
  174.         }
  175.         else    {
  176.         register double    *re_c = re_i + f*vsize,
  177.                 *im_c = im_i + f*vsize,
  178.                 *re_out = re_o, *im_out = im_o;
  179.             for (i=0; i<row; i++){
  180.                 for (j=0; j<dimen1len; j++){
  181.                 *re_out++ = *re_c++; /* if nothing re_c */
  182.                 *im_out++ = *im_c++; /* will be 100 times fast */
  183.                 }
  184.                 re_out += hcols-1;
  185.                 im_out += hcols-1;
  186.             }
  187.         }
  188.         dvrft_2d(re_o, im_o, logrows, logcols);
  189.  
  190.         if (uimg.o_form == IFMT_BYTE){
  191.             dtob(re_o, im_o, fsize, 1./((double)fsize*(uimg.in_form&1?frm:1)));
  192.             temp = im_o;
  193.         }
  194.         else
  195.         {
  196.         register double    total = 1. / (fsize*(uimg.in_form&1?frm:1));
  197.         register float    *op = (float*)(temp = re_o);
  198.             for (i=0; i<fsize; i++)
  199.             op[i] = temp[i] * total;
  200.         }
  201.         fwrite(temp, uimg.pxl_out, fsize, out_fp);
  202.     }
  203. }
  204. else    {
  205. float    *temp, *re_i, *im_i, *cvt,
  206.     *re_o = (float*)nzalloc(fsize, sizeof(float)<<1, "ibuf"),
  207.     *im_o = re_o + fsize;
  208.  
  209.     if (!dimens){
  210.         re_i = (float*)nzalloc(frm*vsize, sizeof(float)<<1, "obuf");
  211.         im_i = re_i + frm*vsize;
  212.         cvt = re_o;
  213.     }
  214.     else    cvt = (float*)nzalloc(fsize, sizeof(*cvt)<<1, "cvt");
  215.  
  216.     if (!dimens){
  217.         w_load(1<<logfrms);
  218.         if (vflag){
  219.             for (f=0; f<frm; f++){
  220.             upread(re_i+f*vsize, sizeof(*re_i), vsize, stdin);
  221.             upread(im_i+f*vsize, sizeof(*im_i), vsize, stdin);
  222.             }
  223.             for (i=frm*vsize; i--;)
  224.             im_i[i] = -im_i[i];
  225.         }
  226.         else    {
  227.         register float    *fop,
  228.             *re_c = re_i,
  229.             *im_c = im_i;
  230.             for (f=0; f<frm; f++){
  231.                 fop = cvt;
  232.                 i = upread(fop, sizeof(*fop)<<1, vsize, stdin);
  233.                 if (i != vsize)
  234.                 syserr("f%d [%d] read %d", f, vsize<<3, i);
  235.                 for (i=vsize; i--;){
  236.                 *re_c++ = *fop++;
  237.                 *im_c++ = -*fop++;
  238.                 }
  239.             }
  240.         }
  241.         for (i=0; i<vsize; i++)
  242.             vfftn(re_i+i, im_i+i, logfrms, vsize);
  243.     }
  244.  
  245.     for (f=0; f<frm; f++){
  246.         if (dimens){
  247.         register float    *re_out = re_o, *im_out = im_o;
  248.             temp = cvt;
  249.             upread(temp, sizeof(*temp)<<1, vsize, stdin);
  250.             for (i=0; i<row; i++){
  251.                 for (j=0; j<dimen1len; j++){
  252.                 *re_out++ = *temp++;
  253.                 *im_out++ = -*temp++;
  254.                 }
  255.                 re_out += cln - dimen1len;
  256.                 im_out += cln - dimen1len;
  257.             }
  258.         }
  259.         else    {
  260.         register float    *re_c = re_i + f*vsize,
  261.                 *im_c = im_i + f*vsize,
  262.                 *re_out = re_o, *im_out = im_o;
  263.             for (i=0; i<row; i++){
  264.                 for (j=0; j<dimen1len; j++){
  265.                 *re_out++ = *re_c++; /* if nothing re_c */
  266.                 *im_out++ = *im_c++; /* will be 100 times fast */
  267.                 }
  268.                 re_out += hcols-1;
  269.                 im_out += hcols-1;
  270.             }
  271.         }
  272.         vrft_2d(re_o, im_o, logrows, logcols);
  273.  
  274.         if (uimg.o_form == IFMT_BYTE){
  275.             ftob(re_o, im_o, fsize, 1./(fsize*(uimg.in_form&1?frm:1)));
  276.             temp = im_o;
  277.         }
  278.         else
  279.         {
  280.         register float    total = 1. / (fsize*(uimg.in_form&1?frm:1));
  281.         temp = re_o;
  282.             for (i=0; i<fsize; i++)
  283.             temp[i] *= total;
  284.         }
  285.         fwrite(temp, uimg.pxl_out, fsize, out_fp);
  286.     }
  287. }
  288. }
  289.  
  290. ftob(ibp, obp, n, scale)
  291. register float*    ibp;
  292. register byte    *obp;
  293. register int    n;
  294. register float    scale;
  295. {
  296. #ifdef    _DEBUG_
  297. message("%d FtoB, scale=%f\n", n, scale);
  298. #endif
  299. if (cut) while (n--)
  300.     *obp++ = *ibp++ * scale;
  301. else for (;n--; obp++)    {
  302.     register int    tmp = *ibp++ * scale;
  303.     if (tmp < 0)
  304.         *obp = 0;
  305.     else if (tmp>255)
  306.         *obp = 255;
  307.     else    *obp = tmp;
  308. }
  309. }
  310.  
  311. dtob(ibp, obp, n, scale)
  312. register double    *ibp;
  313. register byte    *obp;
  314. register int    n;
  315. register double    scale;
  316. {
  317. register int    tmp;
  318.  
  319. #ifdef    _DEBUG_
  320. message("%d DtoB, scale=%lf\n", n, scale);
  321. #endif
  322. for (; n--; obp++)    {
  323.     tmp = *ibp++ * scale;
  324.     if (tmp < 0)
  325.         *obp = 0;
  326.     else if (tmp>255)
  327.         *obp = 255;
  328.     else    *obp = tmp;
  329. }
  330. }
  331.  
  332. dtof(ibp, obp, n, scale)    /* for non power of 2 processing */
  333. register double    *ibp;
  334. register float    *obp;
  335. register int    n;
  336. register double    scale;
  337. {
  338. #ifdef    _DEBUG_
  339. message("%d DtoF, scale=%lf\n", n, scale);
  340. #endif
  341. while (n--)
  342.     *obp++ = *ibp++ * scale;
  343. }
  344.